VM0048-compliant delineation of valid leakage area for the Gola REDD+ Forest Carbon Project project area expansion
country = sf::read_sf("./data/AOI/liberia_boundary_national.shp") |>sf::st_transform(32629)
counties = sf::st_read("./data/AOI/places_poly_county.shp") |>sf::st_transform(32629)
Reading layer `places_poly_county' from data source
`/Users/seamus/repos/rspb-redd-risk-new/data/AOI/places_poly_county.shp'
using driver `ESRI Shapefile'
Simple feature collection with 16 features and 3 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -11.50675 ymin: 4.353908 xmax: -7.367323 ymax: 8.551925
Geodetic CRS: WGS 84
jurisdiction = counties |>dplyr::filter(name=="Grand Cape Mount County"|name=="Gharpolu County")
jurisdiction$name = 'Grand Cape Mount & Gharpolu Counties'
aoi = sf::read_sf("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/AOI/Archive/ProjectArea_CF-Expansion_051525/updated_areas.shp") |>
sf::st_make_valid() |>
sf::st_transform("EPSG:32629") |>
sf::st_cast("MULTIPOLYGON") |> sf::st_as_sf()# |> dplyr::select("Name")
aoi$area_ha = round(as.numeric(sf::st_area(aoi) * 0.0001, 4))
aoi |> sf::st_drop_geometry() |> janitor::adorn_totals() |>
flextable::flextable() |>
flextable::fontsize(size=8,part="all") |>
flextable::autofit()
Name | Shape_Leng | Shape_Area | area_ha |
|---|---|---|---|
Lukasu | 138,549.68 | 469,973,056 | 46,997 |
Mbarma CF | 91,209.87 | 443,230,136 | 44,323 |
Lower Sokpo CF | 59,692.55 | 140,691,204 | 14,069 |
Upper Sokpo CF | 62,831.33 | 109,225,715 | 10,923 |
Zuie CF | 92,221.99 | 362,360,113 | 36,236 |
Gola Forest National Park | 182,951.78 | 890,423,315 | 89,042 |
Tonglay CF | 76,824.36 | 195,953,085 | 19,595 |
Norman CF | 77,870.76 | 123,489,968 | 12,349 |
Foya | 351,458.82 | 1,048,825,566 | 104,885 |
Kposo | 77,241.08 | 307,632,539 | 30,763 |
Gola Konneh CF | 117,719.73 | 531,018,504 | 53,102 |
Gola Konneh CF (ext. S) | 56,323.80 | 108,968,871 | 10,897 |
Gola Konneh CF (ext. N) | 96,399.13 | 285,898,112 | 28,590 |
Lower Sokpo CF (ext.) | 45,714.33 | 47,712,283 | 4,771 |
Maima CF (ext.) | 59,344.81 | 139,163,854 | 13,916 |
Zuie CF (ext.) | 17,920.40 | 14,953,753 | 1,495 |
Yangaya | 46,462.77 | 58,336,055 | 5,834 |
Koninga | 71,795.27 | 317,410,805 | 31,741 |
Lobarsu | 62,061.63 | 165,254,549 | 16,548 |
Bade | 99,147.36 | 513,666,252 | 51,367 |
No-mans land | 209,804.70 | 394,466,117 | 39,447 |
Lower Guma | 56,222.96 | 142,192,367 | 14,219 |
Central Guma | 31,669.60 | 38,501,602 | 3,850 |
Upper Guma | 54,605.75 | 137,757,843 | 13,771 |
Hassala | 62,847.92 | 145,701,577 | 14,570 |
Buluyeama | 47,779.35 | 100,212,025 | 10,021 |
Hembeh | 199,395.60 | 461,835,431 | 46,184 |
Total | 2,546,067.34 | 7,694,854,697 | 769,505 |
aoi_union = sf::st_transform(aoi, 32629) |> sf::st_union() |> sf::st_make_valid()
leakage_buffer = sf::st_buffer(aoi_union, dist = 5500, endCapStyle="ROUND") |>
sf::st_as_sf()
leakage_buffer = concaveman::concaveman(leakage_buffer, concavity=5) |>
sf::st_zm() |> sf::st_make_valid()
leakage_belt_whole = sf::st_buffer(leakage_buffer, dist = 4500, endCapStyle="ROUND") |>
sf::st_make_valid()
tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt_whole) +
tmap::tm_polygons(col="orange",fill="orange",fill_alpha=0.3,lwd=1)+
#tmap::tm_add_legend(type="lines",col="orange",labels="Leakage Belt (10km)") +
tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1, col="white") +
tmap::tm_basemap("Esri.WorldImagery")
Reading layer `LeakageBelt_10k-Radius_UnFiltered-Whole' from data source `/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -11.31102 ymin: 7.024169 xmax: -9.847199 ymax: 8.24916
Geodetic CRS: WGS 84
leakage_belt = sf::st_difference(leakage_belt_whole, st_union(st_combine(aoi_union)))
leakage_belt$area_ha = round(as.numeric(sf::st_area(leakage_belt) * 0.0001, 4))
leakage_belt = sf::st_intersection(country, leakage_belt)
tmap::tmap_mode("view")
tmap::tm_shape(leakage_belt) +
tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.3)+
#tmap::tm_add_legend(type="lines",col="yellow",labels="Leakage Belt (10km)") +
tmap::tm_shape(aoi) + tmap::tm_borders(lwd=0.4, col="white") +
tmap::tm_text(text="Name", size=0.5, col="white") +
tmap::tm_basemap("Esri.WorldImagery")
# save locally
#sf::wt_write(leakage_belt, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/
# Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered.zip")
#sf::st_write(leakage_belt_whole, "OneDrive.../20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered-Whole.shp")
Reading layer `leakage_belt_10km_liberia_ext' from data source
`/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Archive/LeakageBelt_10k-Radius_UnFiltered/leakage_belt_10km_liberia_ext.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N
Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 245117.9 ymin: 776808.3 xmax: 406590.1 ymax: 908414.4
Projected CRS: WGS 84 / UTM zone 29N
Name Shape_Leng Shape_Area
1 Leakage Belt 10km Radius 138549.7 469973056
geometry area_ha
1 MULTIPOLYGON (((245306.7 79... 298838
roads_one = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp")
roads_two = sf::st_read("~/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp")
# we have simplify mask shapefiles and split them up to shorten computing
# time & avoid crashing. See option for "harsh" simplification on line 163
roads_one_simplified = roads_one |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |>
rmapshaper::ms_simplify(keep=0.5)
roads_two_simplified = roads_two |> sf::st_make_valid() |> sf::st_cast("MULTILINESTRING") |>
rmapshaper::ms_simplify(keep=0.5)
# bigger file needs more simplificaiotn
roads_one_simplified_harsh = rmapshaper::ms_simplify(
roads_one_simplified, keep=0.01)
# now apply buffer operation, but note this takes time. Its
# advised processing inputs as much as possible before running
roads_one_buffer = sf::st_buffer(
roads_one_simplified_harsh,
dist = 10000,
nQuadSegs = 5,
endCapStyle="ROUND",
joinStyle = "ROUND",
mitreLimit = 1,
singleSide = FALSE
)
roads_two_buffer = sf::st_buffer(
roads_two_simplified,
dist = 10000,
nQuadSegs = 5,
endCapStyle="ROUND",
joinStyle = "ROUND",
mitreLimit = 1,
singleSide = FALSE
)
# Combine, dissolve and cast to single feature
roads_mask = sf::st_combine(roads_one_buffer, roads_two_buffer) |>
sf::st_union() |> sf::st_cast("POLYGON")
# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
tmap::tm_shape(roads_one_simplified_harsh) + tmap::tm_lines(lwd=2, col="red") +
tmap::tm_shape(roads_two_simplified) + tmap::tm_lines(lwd=2, col="green") +
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=1, col="pink") +
#tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_scale_bar(position = c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))
# Save output to MASKS folder and purge memory
#sf::st_write(roads_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASKS/LeakageMask_Roads_10km-Buffer_051625.shp", delete_dsn=T)
Reading layer `roads_simplified_one' from data source
`/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_one.shp'
using driver `ESRI Shapefile'
Simple feature collection with 842 features and 16 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 237648.4 ymin: 761962.5 xmax: 400449.6 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N
Reading layer `roads_simplified_two' from data source
`/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/ROADS/Archive/roads_simplified_two.shp'
using driver `ESRI Shapefile'
Simple feature collection with 255 features and 16 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 237648.4 ymin: 762265 xmax: 400334.4 ymax: 923370.5
Projected CRS: WGS 84 / UTM zone 29N
Reading layer `Road_Mask_10km-Proximity_051625' from data source
`/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Roads_10km-Buffer_051625/Road_Mask_10km-Proximity_051625.shp'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 17 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 227660.4 ymin: 752279.5 xmax: 410286 ymax: 933257.1
Projected CRS: WGS 84 / UTM zone 29N
# import inputs
wetlands = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Wetlands/GLWD_EPSG32629.tif")
protected_areas = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp")
Reading layer `WDPA_Mar2025_Public_32629_GOLA' from data source
`/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/HABITAT/Protected Areas/Archive/WDPA_Mar2025_Public_32629_GOLA.shp'
using driver `ESRI Shapefile'
Simple feature collection with 6 features and 30 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 240189.7 ymin: 751293.6 xmax: 411102.7 ymax: 895314.1
Projected CRS: WGS 84 / UTM zone 29N
leakage_belt_crop = sf::st_transform(sf::st_as_sf(leakage_belt_whole), 32629) |> terra::vect()
wetlands = terra::crop(wetlands, leakage_belt_crop, mask=T)
# tidy labeling
code_dict_2 <- data.frame(
id = c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31),
label = c(
"Freshwater lake", # 1
"Large river", # 4
"Small streams", # 7
"Riverine, regularly flooded, forested", # 10
"Riverine, seasonally flooded, forested", # 12
"Riverine, seasonally saturated, forested", # 14
"Riverine, seasonally saturated, non-forested", # 15
"Palustrine, seasonally saturated, forested", # 18
"Ephemeral, forested", # 20
"Ephemeral, non-forested", # 21
"Tropical peatland, forested", # 26
"Other coastal wetland" # 31
))
levels(wetlands) <- code_dict_2
wetlands[wetlands == 0] <- NA
# derive wetland mask
wetlands_mask <- wetlands
wetland_classes <- c(1, 4, 7, 10, 12, 14, 15, 18, 20, 21, 26, 31)
terra::values(wetlands_mask) <- ifelse(terra::values(wetlands) %in% wetland_classes, 1, NA)
# save locally for faster computing
#raster::writeRaster(wetlands_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Wetland-GLWD_051625.tif", overwrite=T)
#sf::st_write(protected_areas, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_ProtectedAreas_WDPA_051625.shp", delete_dsn=T)
# skipping these operations here (est. time 12 mins)
DEM = terra::rast("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/DEM/DEM_EPSG32629.tif")
# derive slope percentage from degree
slope_degrees = terra::terrain(DEM, v="slope", unit="degrees")
slope_percent = tan(slope_degrees * (pi / 180)) * 100
slope_percent = terra::clamp(slope_percent, 0, 100)
slope_invalid = slope_percent > 10
slope_invalid[slope_invalid == 0] <- NA
slope_mask = terra::as.polygons(slope_invalid, dissolve=T)|>sf::st_as_sf()|>sf::st_union()
# save locally & reload
sf::st_write(slope_mask, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/LeakageMask_Slope10%-Invalid_051625.zip", delete_dsn=T)
slope_mask = sf::st_read("/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/MASK/Archive/LeakageMask_Slope10%-Invalid_051625/slope_poly_simplified.shp")
# Visual check
tmap::tmap_mode("view")
tmap::tm_shape(roads_mask) + tmap::tm_borders(lwd=0) +
tmap::tm_shape(wetlands) + tmap::tm_raster(col.legend = tm_legend("Wetlands (GLWD")) +
tmap::tm_shape(leakage_belt) + tmap::tm_polygons(col="yellow",fill="yellow",fill_alpha=0.4, lwd=1.5)+
tmap::tm_shape(protected_areas) + tmap::tm_polygons(fill="ORIG_NAME", fill.legend = tm_legend("Protected Areas (WDPA)")) +
tmap::tm_shape(aoi) + tmap::tm_borders(lwd=1.5, col="red") +
tmap::tm_text(text="Name", size=0.3, col="black") +
tmap::tm_shape(roads_mask) + tmap::tm_borders(col="purple") +
#tmap::tm_graticules(lines=T,labels.rot=c(0,90),lwd=0.2) +
tmap::tm_scalebar(position = c("RIGHT", "BOTTOM"), text.size = .5) +
tmap::tm_compass(color.dark = "gray60", text.color = "gray60", position = c("left", "top"))
# clip
roads_leakage_mask = sf::st_intersection(leakage_belt, roads_mask)
slope_leakage_mask = sf::st_intersection(leakage_belt, slope_mask)
wetlands_leakage_mask = sf::st_intersection(leakage_belt, wetlands_mask)
protected_areas_leakage_mask = sf::st_intersection(leakage_belt, protected_areas)
# merge. Btw actually makes better sense to keep these seperate.
# They are easier to operate seperately due to their size and linework.
leakage_mask_a = sf::st_union(roads_leakage, slope_leakage_mask)
leakage_mask_b = sf::st_union(leakage_mask_a, wetlands_leakage_mask)
leakage_mask_c = sf::st_union(leakage_mask_b, protected_areas_leakage_mask)
# save
sf::st_write(leakage_mask_c, "/Users/seamus/Library/CloudStorage/OneDrive-WinrockInternationalInstituteforAgriculturalDevelopment/20087 - RSPB Gola Feasibility/Deliverables/Spatial Data/LEAKAGE/Leakage Masks/", delete_dsn=T)
road_count_whole = sf::st_intersection(road, leakage_belt_whole)
road_count = sf::st_intersection(road, leakage_belt)
road_length_whole = sum(sf::st_length(road_count_whole)) + sum(sf::st_length(road_count_whole))
road_length = sum(sf::st_length(road_count)) + sum(sf::st_length(road_count))
road_length_whole
road_length
waterways_count_whole = sf::st_intersection(waterways, leakage_belt_whole)
waterways_count = sf::st_intersection(waterways, leakage_belt)
waterways_length_whole = sum(sf::st_length(waterways_count_whole))
waterways_length = sum(sf::st_length(waterways_count))
waterways_length_whole
waterways_length
places_count_whole = sf::st_intersection(places, leakage_belt_whole)
places_count = sf::st_intersection(places, leakage_belt)
places_count_whole
places_count